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The Berezinskii-Kosterlitz-Thouless (BKT) phase transition drives the unbinding of topological 
defects in many two-dimensional systems. In the two-dimensional Coulomb gas, it corresponds to an 
insulator-conductor transition driven by charge deconfinement. We investigate the global topological 
properties of this transition, both analytically and by numerical simulation, using a lattice-field 
description of the two-dimensional Coulomb gas on a torus. The BKT transition is shown to be an 
ergodicity breaking between the topological sectors of the electric field, which implies a definition of 
topological order in terms of broken ergodicity. The breakdown of local topological order at the BKT 
transition leads to the excitation of global topological defects in the electric field, corresponding to 
different topological sectors. The quantized nature of these classical excitations, and their strict 
suppression by ergodicity breaking in the low-temperature phase, afford striking global signatures 
of topological-sector fluctuations at the BKT transition. We discuss how these signatures could be 
detected in experiments on, for example, magnetic films and cold-atom systems. 


I. INTRODUCTION 

Topological physics [T] emerges in many condensed- 
matter systems, including superfluids and super¬ 
conductors [SHU, topological insulators [5], exciton- 
polariton condensates |^, and magnetic textures such 
as skyrmions [7]. Among two-dimensional systems, a 
prototypical application of topology concerns the quan¬ 
tum Hall effect in the two-dimensional electron gas i- 
OTO], while many other examples relate to the physics of 
topological defects identified by Berezinskii, Kosterlitz 
and Thouless (BKT) [TTHI^ . These include Josephson 
junction arrays [IIHIT], films composed of Bose-Einstein 
condensates [nmn], superfluid films |20| . liquid-crystal 
and polymer films [U], superinsulators [1HI23], and mag¬ 
netic films and layers [24H27] . In such systems, the BKT 
phase transition drives the thermal dissociation of bound 
pairs of local topological-defects [TTHISI HH] • The idea of 
a topological defect (defined in the footnote [5^) is in¬ 
deed one of the most basic and important applications of 
topology in condensed-matter physics [50] . 

An important discovery of BKT and later authors nu¬ 
lla ng is that the defect-mediated transition of the plane 
rotator (or 2D-XY) model and its analogues can be 
mapped to the insulator-conductor transition of a two- 
dimensional Coulomb gas m- The long-range Coulomb 
interactions emerge from a purely local Hamiltonian, so 
that the mapping at the microscopic level, although com¬ 
plete [Sg, is far from transparent. However, as Maggs 
and co-workers [33H38| have shown in three dimensions, 
a Coulomb fluid can be transformed into a local prob¬ 
lem by using an electric-field representation and intro¬ 
ducing a freely fluctuating auxiliary gauge field. Fol- 
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lowing this work, it is straightforward to show that the 
XY Hamiltonians that admit a BKT transition map on 
to this generalized electrostatic problem in two dimen¬ 
sions. A practical consequence of the phase-space ex¬ 
tension to a fluctuating auxiliary gauge field is the de¬ 
velopment of purely local algorithms for the simulation 
of Coulomb fluids in both three [35U38] and two dimen¬ 
sions [39|, which circumvent the technical difficulties as¬ 
sociated with long-range interactions. In particular, the 
logarithmic potential that governs charge-charge interac¬ 
tions in the two-dimensional Coulomb gas is dealt with 
locally, allowing a new approach to the efficient simula¬ 
tion of two-dimensional Coulombic systems. 

In this paper, we exploit these developments to formu¬ 
late and simulate a lattice-field description of the two- 
dimensional Coulomb gas for the purpose of investigating 
the topological properties of the BKT transition. The 
BKT transition is topological in the sense that it sep¬ 
arates a topologically ordered phase from a disordered 
one. Topological order in this context means that the 
local topological defects (charges in the two-dimensional 
Coulomb gas) are confined. Vallat and Beck [32] con¬ 
sidered the two-dimensional XY model on a torus, and 
showed how a winding field can be associated with the 
global topology of the system. In the high-temperature 
phase, where charge is deconfined, non-zero values of 
this winding field define global topological defects that 
are distinct from the local topological defects driving the 
BKT transition. Here we show that the lattice-field de¬ 
scription naturally lends itself to classifying and investi¬ 
gating this property. In this paper, we treat only the two- 
dimensional Coulomb gas, but in a further publication we 
will extend our analysis to the case of two-dimensional 
XY models on the torus. 

Our key observation is that the topology of the torus on 
which the Coulomb gas is placed generates a multiplicity 
of states in the lattice electric-field representation that 
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are equivalent for charge configurations but not energeti¬ 
cally degenerate. Given an arbitrary charge distribution, 
one is at liberty to add an integer multiple of some con¬ 
stant to each component of the harmonic mode of the 
electric field while leaving the charge distribution un¬ 
changed. This global topology associated with the BKT 
transition describes the winding of charges around the 
torus. In the high-temperature phase, charge deconfine¬ 
ment allows for fluctuations in the winding component 
of the harmonic mode, which can be classified as differ¬ 
ent topological sectors. Below the transition, however, 
the binding of charge pairs causes the winding compo¬ 
nent to be zero. Topological-sector fluctuations in the 
electric field therefore mark the appearance of the high- 
temperature, topologically disordered phase at the BKT 
transition. 


The present study of topological-sector fluctuations in 
the two-dimensional Coulomb gas may be compared to 
previous studies on the three-dimensional Coulomb phase 
of spin-ice materials and models In spin ice, the 

onset of topological-sector fluctuations is shown to sig¬ 
nal a Curie law crossover [40] for the zero-field suscep¬ 
tibility and a Kasteleyn transition in the presence of an 
applied held mma. Our study of the BKT transition 
reveals aspects of topological-sector huctuations that are 
not found in either of these established cases. For exam¬ 
ple, by our analysis, the two-dimensional Coulomb gas 
should be considered to present an ergodicity-breaking 
transition to a topologically ordered phase in the absence 
of an applied held, whereas spin ice has no equivalent 
phase. 

The paper is organized as follows. In Section [TT| 
we introduce the lattice-held representation of the two- 
dimensional Coulomb gas on a torus and use this to dehne 
the partition function and the topological sectors of the 
electric held. We use numerical simulations to demon¬ 
strate that topological-sector huctuations appear in the 
high-temperature (conducting) phase but n ot in the low- 
temperature (insulating) phase. In Section III we show 
that the reason for the strict suppression of topological- 
sector huctuations in the low-temperature phase is er- 
godicity breaking at the transition. A hnite-size scaling 
analysis, given in Section EYI conhrms that in the ther¬ 
modynamic limit, ergodicity is broken precisely at Tbkt- 
Conclusions and comparisons with experimental systems 
are discussed in Section M 


II. TOPOLOGICAL-SECTOR FLUCTUATIONS 


Using the unit system dehned in Appendix [A] we for¬ 
mulate the two-dimensional Coulomb gas using discrete 
vector calculus on a square lattice with periodic bound¬ 
ary conditions (PBCs) applied. The PBCs enforce the 
toroidal topology but not the curvature of a true torus. 
All functions are dehned to be the discrete counterparts 
of smooth vector helds jUj, and any lattice vector held 


F is dehned component-wise via [H] 

F(x) := F, (x -k -k F), (x -k By, (1) 

where x is any lattice site and is the unit vector 
in the x/y direction. The operators V and V are the 
forwards and backwards hnite-difference operators on a 
lattice, resi^ectiyely, and the lattice Laplacian is dehned 
by := V • V jUj. The most general electric held E 
may be Helmholtz decomposed into the sum of a Poisson 
(divergence-full) component — Vi)), a rotational compo¬ 
nent E and a harmonic component E: 

E(x) = -V<))(x)-kE(x)-kE. (2) 

This electric held is the most general solution to Gauss’ 
law on a lattice: 


V • E(x) = p(x)/eo. (3) 

Here, p(x) := qm(x)/a^ is the charge density at each 
lattice site x, q is the elementary charge, the integer m 
denotes the charge species, a is the lattice spacing and cq 
is the electric permittivity of free space (see Appendix]^ . 
Using the held E adds an auxiliary held E to the usual 
solution of electrostatics, as in the electrostatic model 
of Maggs and Rossetto (MR) [^. This allows us to 
simulate the physics of Coulombic interactions on a lat¬ 
tice via local electric-held updates, avoiding the need to 
treat computationally intensive long-range interactions, 
as outlined in Appendix |Bj 

The validity of introducing the auxiliary held is seen in 
the context of the separability of the partition function 
into its Coulombic and auxiliary components: the auxil¬ 
iary held contributes to the internal energy of the system, 
but it is statistically independent of the Coulombic ele¬ 
ment. In Appendix we give a full description of the 
algorithm and a derivation of the partition function for 
the Coulomb gas of multi-valued charges. 

The internal energy of the electric helds corresponding 
to a given charge and auxiliary-held conhguration is given 
by 


C^o = ^E|EWI'- (4) 

where D is the set of all lattice points. To represent the 
Coulomb gas in the grand canonical ensemble, we add a 
core-energy term Ucoie, given by 

Ucore ^ E 
xeD 

where ec(m) is the core-energy constant of each charge 
mq, and edm) = ed—m) since charges are excited to 
the vacuum in neutral pairs. The grand-canonical en¬ 
ergy of the system U = Uq + Ucore may be expanded by 
combining Eqs. (i), 0 and 0 to give a sum of terms 
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arising from the different field components, which add to 
the core energy: 


U = Uself + Ulnt + C^Rot + C^Harm + Ucore- ( 6 ) 

Here, respectively, HRot := X^xs-D and 

C^Harm := eo-^^|Ep/2 are the auxiliary-field and 

harmonic-mode components of the grand-canonical en¬ 
ergy, and Hseif and Hint are the self-energy and Coulom- 
bic charge-charge interaction components. As outlined 
in detail in Appendix]^ the latter two components may 
be expressed in terms of the lattice Green’s function G, 
according to Hseif := a'‘ G(0) Gint := 

«'‘Ex.5^x,GDf^(^*)G'(xi,Xj)p(xj)/2eo, where G(0) := 
G(x, x). Note that, while Hint can be negative, the sum 
Gseif + Gint is necessarily > 0, as it arises from the term 
in |V(/)p. 

Using the above results, we may define the chemical 
potential for the introduction of a charge mq: 


■ — 


G(0) 

Co 


ec{m) 


2 9 
2 


( 7 ) 


In the following, we specialize to a Coulomb gas of n pairs 
of elementary charges of chemical potential /i := fii, by 
setting ec{m = 0, ±1) = 0 and €c{m ^ 0,±1) = oo [45]. 

The harmonic mode E is a uniform field found by av¬ 
eraging the total electric field E(x) over x. In a simply 
connected space, the average field may be conveniently 
related to the average polarization P arising from an ef¬ 
fective surface charge distribution by E = — P/eq. For 
a charge-neutral system in a simply connected space, 
P := is invariant with respect to the ori¬ 

gin shift X I—> x-I-xq, and is therefore origin-independent. 
The situation is more complicated on a toroidal surface 
as E can also depend on a harmonic-field component that 
corresponds to a charge winding around the torus, and it 
is necessary to adopt a convention to define distances be¬ 
tween points (the concepts ‘close together’ and ‘far apart’ 
are ambiguous on a torus). In Appendix]^ we show in 
detail how it is possible to define origin-independent po¬ 
larization Ep and winding E^ components of the har¬ 
monic mode such that 


E — Ep + Ew. (8) 

Here, 

Ew = -^w, (9) 

Leo 


wind around the torus in opposing directions before as¬ 
suming its original configuration. When a single charge 
winds around the torus in the x/y direction, the x/y 
component of the harmonic mode of the electric field 
increases by Eg/Leo- As shown in Appendix 
the lowest-energy harmonic mode that describes an ar¬ 
bitrary charge distribution is therefore an element of the 
set {—q/2Leo, q/2Leo] and is defined as the polarization 
component in Eq. (10) by applying modular arithmetic 


to E. The remainder is the winding component. The 
modulo operation removes any need for a ‘distances’ con¬ 
vention to define the polarization component, as well as 
any origin dependence of the field components (see Ap¬ 
pendix]^ for further details). 

With these results we may use the integer-valued wind¬ 
ing field w to define the topological sector of the system 
as the number of times charges wind around the torus 
in the x and y directions, with all non-trivial topological 
sectors given by w ^ 0. The topological sector of the sys¬ 
tem changes any time a charge pair unbinds and winds 
around the torus and hence thermal fluctuations of the 
topological sector are closely related to the unbinding of 
charge pairs, as elucidated further below. 

The statistical mechanics of the topological-sector fluc¬ 
tuations may now be formulated by considering how the 
polarization and winding components of the harmonic 
mode enter the lattice partition function. As shown in 
Appendix B, the partition function splits into two sta¬ 
tistically independent components. One component is 
the Coulombic partition function .^coui and contains all 
information about the charge-charge correlations, while 
the other is the auxiliary-field partition function .Zroi 
and contains all information about the auxiliary field: 
the auxiliary field can freely fluctuate without affecting 
charge-charge correlations (see Appendix [b|) . The parti¬ 
tion function is written as 


E — .^Coul-^Rot; 

where Zcoui is given by 


( 11 ) 


Zconl = ^ exp 
{p(x)} 

X exp 


2eo 


Y P(x*)G(x„Xj)p(xj) 


i2/3eO|^ q ,2 



( 12 ) 


where w is an integer-valued winding field chosen such 
that 


E, 


'P,xly 


g g 
2Heo 2Leo 


( 10 ) 


and L is the lattice length. 

This decomposition of the harmonic field has the fol¬ 
lowing interpretation. A charge pair may unbind and 


Here, /3 := l/k^T is the inverse temperature and the sum 

^{p(x)} ■~^{a‘^p(x)G{0,±(3}}- 

The first exponential of Eq. (121 describes the anhar- 
monic charge-charge interactions, the second describes 
the polarization and winding state of the system, and 
the third describes the sum of the self-energies associated 
with each elementary charge. The sum over the winding 
field w is necessitated by the degeneracy of the harmonic 
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Monte Carlo time steps /10^ 


FIG. 1. Topological-sector fluctuations of lattice electric fields 
in the two-dimensional Coulomb gas on a torus. Shown is the 
a;-component of the normalized total harmonic mode LE^I^-e 
and winding field L.Bw,a:/27r versus Monte Carlo time for an 
L X L system of linear size Z/ = 16atT = 1.34 (top) and 
T = 2.0 (bottom). The system was simulated using the MR 
algorithm with local moves only. At the lower temperature 
(top) harmonic-mode fluctuations are finite (black) but there 
is no winding-field component (blue), while at the higher tem¬ 
perature (bottom) the winding-field component becomes fi¬ 
nite, indicating topological-sector fluctuations. 


mode of the electric field: infinitely many topological sec¬ 
tors describe any given charge configuration. In addition 
to the local updates of the MR algorithm (see Appendix 
0. we also consider global updates, which correspond to 
independently sampling this winding field. 

Henceforth, we set the elementary charge q = 27r, the 
lattice spacing a = 1, the electric permittivity eo = 1, and 
Boltzmann’s constant /cb = 1. The choice q = 27r recog¬ 
nizes the standard BKT theory, where a charge emerges 
as a local 27r winding in an associated lattice field, such 
as the spin differences in the 2D-XY model [55]. 

The BKT transition drives the deconfinement of charge 
pairs in the two-dimensional lattice Coulomb gas, which 
generates topological-sector fluctuations. The transition 


occurs at Tbkt = 1.351 (to four significant figures) [15] 
in the thermodynamic limit [a value specific to a gas 
of elementary charges with edm = 1) = 0], which is 
scaled to higher temperatures in finite-size systems (see 
below). Fig. shows the evolution of the (normalized) 
cc-component of the harmonic mode of a system of linear 
size L = 16, simulated using local moves only (numeri¬ 
cal simulation details are described in Appendix No 
topological-sector fluctuations are visible just below the 
BKT transition temperature Tbkt = 1.351, but they be¬ 
come important at temperatures above Tbkt- 


III. ERGODICITY BREAKING 


A convenient measure of the topological-sector fluctu¬ 
ations is the winding-field susceptibility Xw, defined by 

Xw(T,T) := /3eoT' ((E^) - (Ew)') - (13) 


In a fully ergodic system, Xw is nonzero, even in the ab¬ 
sence of charge fluctuations, as can be seen by limiting 
the Gibbs ensemble contributing to ^coui to field config¬ 
urations of zero charge. In this case it is straightforward 
to show, using Eqs. ([^ and (121, that the constrained 
susceptibility is given by 


^global(j,) ^ ^^^^2 


exp {-Pq'^/2eo) /til? E ... 
1 -I- 4 exp (-/3g2/2eo) + ... 


eo 


exp (-/3gV2eo) , 


(14) 


for k^T q^l2eo. The system-size dependence falls out 
of this expression so that a fully ergodic system would 
show small but finite topological-sector fluctuations in 
the low-temperature phase. 

Assuming local charge dynamics, a topological-sector 
fluctuation requires the separation of a pair of charges 
over a distance greater than T/2 in either the x or the 
y direction [see Eq. ^ and the subsequent discussion]. 
As the charge concentration falls to zero at low temper¬ 
ature, screening becomes negligible and the energy bar¬ 
rier against such configurations diverges logarithmically 
with the linear system size L dniiaini]. As the charge 
concentration increases with temperature, however, en¬ 
tropy and charge screening break down the free-energy 
barrier, making it finite at the BKT transition. Above 
the transition, charge pairs are free to unbind and trace 
closed paths around the torus, giving finite-valued wind¬ 
ing fields, as observed in Eig. [^ In contrast, in the low- 
temperature phase, the probability of separation through 
a distance T/2 becomes strictly zero in the thermody¬ 
namic limit. 

The BKT transition is therefore an ergodicity break¬ 
ing: a change in the phase space explored by a system 
with local dynamics. In detail, it is an ergodicity break¬ 
ing between topological sectors, signalled by the strict 
suppression of topological-sector fluctuations in the elec¬ 
tric field at T < Tbkt- If the dynamics were non-local 































5 



FIG. 2. The susceptibility quotient Xw'^'^VXw* versus temper¬ 
ature for an L X L Coulomb gas of linear size L = 64. In 
the region T < 1.2, the quotient is zero, while for T > 1.6, 
the quotient approaches unity. This divergence between the 
results of the local-update and the all-updates simulations, 
accompanied by striking fluctuations in the intermediate re¬ 
gion, signals an ergodicity breaking as the system is cooled 
through the BKT transition. The line is a guide to the eye. 


(including global updates of the winding component of 
the harmonic mode [55h38) l. Xw would remain finite at 
all temperatures. 

To explore this ergodicity breaking, we have simulated 
the two-dimensional Coulomb gas, first with local field 
updates only, and second with both local and global 
field updates [33HSH]. Corresponding to each case, we 
define the winding-held susceptibilities and Xw'j 

respectively. Differences between and Xw* reflect 

the inability of local moves to explore a fully represen¬ 
tative phase space on the time scale of the simulation. 
To quantify this, we introduce the susceptibility quotient 
Xw'^^VXw^ which may be used to analyse the ergodicity 
of the system. 

Fig. [2] clearly shows that ergodicity is broken in the 
vicinity of the BKT transition. For T > 1.6, = Xw'i 

indicating that the free-energy barrier for a topological- 
sector fluctuation via local moves is small. For T < 1.2, 
the quotient is zero, indicating that the energy barrier 
prevents topological-sector huctuations via local charge 
moves. In between these low- and high-temperature re¬ 
gions there are strong huctuations in the quotient be¬ 
cause charge deconhnement via local updates represents 
increasingly rare events, an inevitable precursor to loss 
of ergodicity. In Section]^ this ergodicity breaking is 
shown to occur precisely at Tbkt in the thermodynamic 
limit. 

Our analysis thus leads to a precise dehnition of 
topological order for the two-dimensional Coulomb gas 
through the ergodic freezing of the topological sector to 
its lowest absolute value. Two-dimensional systems with 
C/(l) symmetry are often associated with an absence of 


an ordering held at hnite temperature m- Here we ex¬ 
plicitly show that, in the case of the BKT transition, the 
ordering of a conventional order parameter is replaced 
by topological ordering through an ergodicity breaking 
between the topological sectors. The topological order is 
directly related to the conhnement-deconhnement tran¬ 
sition of the charges, the local topological defects of the 
electric held. This type of ergodicity breaking is dis¬ 
tinct from either the symmetry breaking that character¬ 
izes a standard phase transition, or that due to the rough 
free-energy landscape that develops at a spin-glass tran¬ 
sition [15] . 


IV. FINITE-SIZE SCALING 

In order to explore the approach to the thermodynamic 
limit, the two-dimensional Coulomb gas was simulated by 
the Monte Carlo method as a function of system size, us¬ 
ing the MR algorithm. The global update was employed 
in order to improve the statistics (numerical simulation 
details are described in Appendix [Pj) . 

Fig. [g shows the simulated winding-held susceptibil¬ 
ity Xw as a function of temperature for L x L Coulomb 
gases of linear sizes between L = 8 and L = 64. There is 
a marked increase in the winding-held susceptibility Xw 
as the system passes through the BKT transition tem¬ 
perature Tbkt = 1.351 US] for all system sizes. Sus¬ 
ceptibility curves for successive values of L intersect at 
temperatures above T = 1.8 and below T = 1.5. Between 
these two temperatures, the winding-held susceptibility 
increases for a given temperature as the linear system size 
L increases. These results are consistent with the hnite- 
size scaling of the BKT transition temperature [24|: as 
the system size decreases the effective transition temper¬ 
ature T* (L) increases. 

Below Tbkt, the probability of a charge pair sepa¬ 
rating over a distance greater than L/2 increases with 
decreasing system size. This, combined with the hnite- 
size transition temperature T*{L) also increasing with 
decreasing system size, results in the winding-held sus¬ 
ceptibility curves for successive values of L intersecting 
in the vicinity of Tbkt- The inset in Fig. shows that 
the low-temperature crossover points of the susceptibil¬ 
ity curves are at T = 1.45, T — 1.40 and T = 1.37 (to 
within estimated error). To extrapolate the trend of the 
data shown in Fig. [^to the thermodynamic limit, we de- 
hne the crossover temperature Tcmss{L) to be the lower 
temperature at which Xw(T) = XwiL/2). 

Fig. a shows the crossover temperature Tcross as a 
function of inverse linear system size 1/T, along with 
straight-line hts to the data. In the thermodynamic limit, 
Tcross extrapolates to the value Tcross = 1.351 to within 
the estimated error of the extrapolation, that is, it ex¬ 
trapolates to the BKT transition temperature [46] : 

Tcross(T —)■ oo) = Tbkt- (15) 
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FIG. 3. The winding-field susceptibility Xw as a function 
of temperature for L x L Coulomb gases of linear size L = 
8, 16, 32 and 64 (nsing local and global MR moves). The 
curves intersect at low and high temperature. Inset: An ex¬ 
panded plot of the data in the region of the low-temperature 
intersections (with error bars representing two standard de¬ 
viations). The indicated crossover temperatures are given 
by rcrosa(T = 16) = 1.45, Tctosb{L = 32) = 1.40 and 
Tcrosa(T = 64) = 1.37 (to within estimated error), based on 
a data fit. 
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FIG. 4. The crossover temperature Tcroaa (black data; left- 
hand y axis) and crossover susceptibility (red data; 

right-hand y axis) as functions of inverse linear system size 
l/L, with error bars representing two standard deviations. 
Lines are weighted (with respect to the error bars) linear- 
regression fits to each data set, from which the y-intercept 
(L —>• oo) was calculated. Tctoss{L —>■ oo) = 1.351(2), equal 
to the BKT transition temperature Tbkt |46| . The crossover 
susceptibility Xw'^°°°(i oo) = 5 x 10“^ with estimated error 
of the same order: there is no measurable difference between 
this quantity and the winding-field susceptibility due to global 
updates only at T = 1.351. 


The l/L scaling of Tcross is unusual for the BKT tran¬ 
sition, for which the finite-size BKT transition temper¬ 
ature typically scales as a simple function of the loga¬ 
rithm of L [24l |49]. However, Minnhagen and Kim [50] 
found that a fourth-order cumulant that measures fluc¬ 
tuations of the helicity modulus in the 2D-XY model 
also scales as l/L: as this closely relates to fluctuations 
in the harmonic-mode susceptibility |32j . it seems likely 
that we are observing the same finite-size scaling here. 
The magnitude of the winding-field susceptibility at the 
crossover points similarly extrapolates 

to ~ 5 X 10“"^ in the thermodynamic limit, with an es¬ 
timated error of the same order. This small number is 
not measurably different to the winding-field susceptibil¬ 
ity due to global moves only, which, at Tbkt, evaluates 
to ^ 5 X 10“® for all system sizes [see Eq. (14)]. The 
inference is that topological-sector fluctuations due to lo¬ 
cal moves only turn on precisely at the universal point 
T'cross(L —?► oo) = Tbkt iu the thermodynamic limit. 
This confirms that topological-sector fluctuations signal 
charge deconfinement and the high-temperature phase of 
the BKT transition. 

Given that the topological-sector fluctuations turn 
on at the temperature at which the system experi¬ 
ences the famous universal jump in the helicity modu¬ 
lus [HUHlEllEni, it is interesting to estimate the con¬ 
tribution that topological-sector fluctuations make to the 
universal jump. To do this, we define the harmonic¬ 
mode susceptibility xg polarization susceptib- 

lity Xp by replacing Ew in Eq. (13) with E and Ep, 


respectively. The helicity modulus is then given by 
T = (1 — Xe/2) [32], so that xg makes a jump of or¬ 

der unity at Tbkt- We find that the ratio (xg~Ap)/xg is 
less than 5 x 10“^ for all T < 1.6 for systems of linear size 
L = 8 to 64, showing that the contribution to the uni¬ 
versal jump from topological-sector fluctuations is very 
small. This reflects the near-cancellation of (E^) with 
the coupling term 2(Ep • E*) in the evaluation of xg, re¬ 
flecting strong correlations between the polarization and 
winding fields at the transition. 


V. CONCLUSIONS 

In conclusion, the BKT transition has long been a 
paradigm for the importance of topological defects in 
condensed-matter physics [T|. Vallat and Beck showed 
that XY-type systems on the torus generate global topo¬ 
logical defects at the BKT transition that reflect the 
toroidal topology [32|. Here we have used lattice-field 
simulations to reveal topological-sector fluctuations in 
the electric field of a two-dimensional lattice Coulomb gas 
on a torus. We have shown how these provide a striking 
and sensitive measure of the topological and ergodicity- 
breaking character of the BKT transition, allowing a pre¬ 
cise definition of topological order in terms of this broken 
ergodicity. 

The topological-sector fluctuations at the BKT tran¬ 
sition are very clearly revealed in the lattice electric 
field description of the two-dimensional Coulomb gas. 
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but we expect them to be equally relevant to any sys¬ 
tem that has a BKT transition. In suitable systems, 
the winding-field susceptibility that signals the onset of 
topological-sector fluctuations will contribute to experi¬ 
mentally measurable responses of the system. For exam¬ 
ple, in a cylindrical or toroidal magnetic film with XY 
symmetry, winding-field fluctuations in the Coulomb gas 
representation correspond to measurable spin configura¬ 
tions in the magnetic representation. As we will show 
in future work m, fluctuations of an appropriate topo¬ 
logical sector accompany the destruction of the finite- 
size magnetization of an XY spin system through vortex 
deconfinement. They could therefore be observable in 
ultrathin ferromagnetic metallic hlms [52] or magnetic 
Langmuir-Blodgett films |551l5Tj . 

Another promising system on which to measure these 
topological-sector fluctuations is the one-dimensional 
quantum lattice Bose gas. When the system is placed 
on a ring, its angular momentum is no longer a good 
quantum number. The angular momentum can there¬ 
fore fluctuate quantum mechanically, and the system 
should undergo a dramatic increase in these fluctua¬ 
tions as it passes through the superfluid - Mott in¬ 
sulator quantum phase transition [551 156] . This dra¬ 
matic increase in the fluctuations corresponds to finite¬ 
valued global topological defects in the quantum system, 
and therefore, via the Feynman path-integral mapping, 
to topological-sector fluctuations in the two-dimensional 
classical lattice Coulomb gas on a torus. Murray et al. 
measured the angular momentum of ring-shaped Bose- 
Einstein condensates via the vortex-density profile of the 
system m- Our measure of the BKT transition could 
therefore correspond to equivalent, experimentally mea¬ 
surable topological-sector fluctuations in the cold-atom 
system [58]. 

Finally, it is worth noting that it is natural to asso¬ 
ciate a conducting phase with the excitation of winding 
fields, as may be seen by considering a loop of wire in 
a changing magnetic field. Recalling that the magnetic 
field does no work on a test charge, the induced elec¬ 
tromotive force must arise from a divergence-free electric 
field running round the loop. The curl of this held obeys 
the Maxwell-Faraday law, V x E = dBldt. Hence, in 
three dimensions, electromagnetic induction provides a 
practical method of exciting topological winding helds 
analogous to those discussed here. 
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Appendix A: Dimensional analysis of the 
two-dimensional Coulomb gas 

In the following, [ • ] denotes the dimensions of some 
quantity, L denotes the dimensions of length, d is the 
spatial dimensionality of the system, and Cq is the vacuum 
permittivity in d—dimensional space. Consider Gauss’ 
law for the MR algorithm, 

V • E(x) = /o(x)/eo, (Al) 

and the dimensions of electric charge density, 

[p(x)] = [g]L-^ (A2) 

which generates 

[E(x)] = [g]L(i-'^)[eo]-\ (A3) 

From a consideration of the exponent of the Boltzmann 

factor (with /3 = l/k^T) we hnd a dimensionless group, 

n E |E(x)P (A4) 

xGD 

[eo] [/3] (A5) 

^ [E(x)] =[g]-iL-M/3]-^ (A6) 

Setting the charge to be dimensionless, it follows that 
[eo] = [/3] (A7) 

and 

[E(x)] = [/3]-iL-i (A8) 

in d = 2. 


Appendix B: The MR electrostatic model and the 
partition function 

The MR electrostatic model is a lattice-held model 
from which it is possible to form the lattice partition 
function of electrostatics. To show this, we describe the 
MR algorithm in terms of microscopic variables that rep¬ 
resent the local held updates. A conjugate lattice D' is 
dehned such that each of its sites is at the centre of each 
plaquette of D. Each site in D' is associated with a real¬ 
valued variable tp whose adjustment corresponds to an 
update of the auxiliary held, while each pair of nearest- 
neighbour sites is associated with an integer-valued vari¬ 
able s whose adjustment corresponds to a charge-hop up¬ 
date. Both sets of variables are subject to PBGs. 

Gomponent-wise, we now dehne the held 

a \ i^(x-l-aei) —(/3(x) + gs(x-|-aei, x) 

2®V a ’ 

(Bl) 




and identify 


E(x) ^ - 
eo 


[A6l]j,(x+fe^) \ 

-[A6»]^(x+ fey) J 


(B2) 


updates to create charge pairs out of the vacuum and 
include the possibility of all integer-valued multiples of 
the elementary charge. Combining Eqs. @ and ( |B2| ), 
the partition function in the microscopic-variable repre¬ 
sentation is given by 


The X coordinates in Eq. (Bl) are in H'; the x coordi- 




-1 

nates in Eq. (B2| are in D. 

A charge hop in the positive x/y direction corresponds 

{«} 

^ Vip exp 

2eo 

Y l<P(x)-</5(x')+gs(x,x')|2 

(x.x'j 

to a decrease/increase in the relevant s variable by an 




amount q, as shown in Fig. (where Sij represents the s ^ ( PUcore ); (B3) 

variable between sites i and j of the conjugate lattice). 

where 


4 4 


1 


€>■ 


-►O 



[vp:=Y[ 

■ r9/2 

/ dp{x) 

J 

^ xGD' 

_J-q/2 


(B4) 


and ■= Z]{s(x,x')GZ}- Here, the grand-canonical 

energy of the Coulombic system U = Uq + C/core is used. 
This representation reproduces Gauss’ law: 


FIG. 5. A charge-hop update in the positive x direction: The 
Sij variable (red arrow) has its value decreased by an amount 
q. The value of the electric field flux Eap (black arrow) flowing 
from site a to site /3 then decreases by g/eo, corresponding to 
a charge-hop update. Red circles represent positive charges; 
white circles represent empty charge sites. 

Fig. [^depicts the microscopic-variable representation 
of the auxiliary-field updates, with an alteration of a 
particular Lp variable rotating the field around its sur¬ 
rounding plaquette. In the figure, the (p variables are 
represented by spin-like arrows in order to emphasize the 
rotation of the electric field. 




FIG. 6. An update of the rotational degrees of freedom of the 
electric field: The value of the p variable at the centre of a 
randomly chosen lattice plaquette decreases by an amount A. 
This rotates the electric flux by an amount A/eo around the 
plaquette, leaving Gauss’ law satisfied. Red arrows represent 
ip variables, black arrows represent the electric field, dashed 
red lines represent the conjugate lattice D', the blue arrow 
represents the direction of the field rotation and grey circles 
represent sites of arbitrary charge. 

With the internal energy of the electric fields corre¬ 
sponding to a given charge and auxiliary-field configura¬ 
tion given by Uq in Fq. it is possible to write the 
partition function in the microscopic-variable represen¬ 
tation. For ease of manipulation, we allow charge-hop 


^ Ad(x) • l(x) = Qr, (B5) 

xGar 

where Qy is the charge enclosed within some subset of the 
lattice T U D,dT (Z D' is the boundary enclosing F, and 1 
traces an anticlockwise path along 9r and has dimensions 
of length. This equation results from the p variables 
cancelling and the s variables being integer valued. It 
follows that 


A6»(x) • l(x) =a^ Y ■ E(x) (B6) 

xGOr xGT 

V • E(x) =p(x)/eo, (B7) 

recovering Eq. (|^, as required. 

The constraints imposed upon the electric field (Gauss’ 
law and the form of the harmonic mode, the latter of 
which is derived in detail in Appendix]^ are combined 
with the grand-canonical energy of the system to write 
the partition function in terms of the electric field. We 
define the set X := such that the partition func¬ 

tion is given by 

Z = \3\ Y E /l^Eexp 5] |E(x)|2 

{/ 9 (x)GX} _ xez> 


X exp i-PUcove) “ P(^)/eo) 



where the functional integral 


VF := TT / dFx{x + aex/2) / dFy{x + aey/2) 

^^20 </]R 


(B9) 
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for any vector field F, and |J| is the Jacobian determi¬ 
nant. 

This partition function may be separated into two com¬ 
ponents by defining the new rotational field 

e(x) := E(x) -H V(/>(x) - E. (BIO) 

The partition function is then given by 

^ = -^Coul ^Rot, (Bll) 


where 


^Coui := ^ exp 

{V20(x)gF} 


/3eoa^ 




xG-D 


X E exp -E|LP-gwo| 
X exp (-PUcore) , 


(B12) 


and 


ZRot ■■= |J| J Veexp 


l3eoa^ 


E 


xGD 


n ['5(^-e(x))]<5( E^C 


xGD 


\xGD 


(B13) 


are the Coulombic and auxiliary-field components of the 
partition function, respectively. Here Y := qZ/coa'^ and 
we have used the fact that all coupling terms in the 
grand-canonical energy sum to zero. The MR algorithm 
therefore reproduces Coulombic physics since the sepa¬ 
ration of the auxiliary-field partition function from the 
Coulombic partition function ensures that the charge- 
charge correlations are independent of the auxiliary field. 

The lattice Green’s function G(x, x') between two 
charge-lattice sites x and x' is defined such that 

V2G(x,x') = -dx.x', (B14) 

where the subscript x denotes with respect to which co¬ 
ordinate system the lattice Laplacian is applied. 

We define the k-space lattice Green’s function G, 

Gx^(k):= Ee-*‘^"G'(x,x'), (B15) 

xGD 


and the set J2k^B ■= YliG{x,y} [Efc.GsJ’ ’^^i^re B, := 
/0±— ±2— ••• ±1^ - 11— =^—1-is the set of 

k-space values in the i direction, and W := '/N. Gom- 
bining Eqs. (B14| and (|B15 ), it then follows that 


E = 2Ee*‘' " [2 - cos(fc,a) - cos(fcj,a)] 

X Gx<k). (B16) 


This is solved by 


_ 2 k' 

Gx'(k) = -f-- ^^—^VkT^ 0 

2 [2 —cos(fca;a) —cos(fcya)] 


(B17) 


where the k = 0 part of the lattice Green’s function is set 
to zero since the harmonic component of E is attributed 
to E. It follows that 


G(x,x') 


1 ^ 

2N 2 — cos{kxa) — cos{kya) 


(B18) 


The internal energy of the Poisson component of the 
electric field is given by 

Gpoisson E (B19) 

xGD 

= - ^ E (B20) 

xGD 

= E E 7'(^*)G'(x„Xj)p(xj), (B21) 

x..x,Gr3 

hence, the Coulombic partition function can be written 
as 


Zcovl\ 


E 

{p(x)GJf} 


2eo 


E p{xi)G{x,,Xj)p{Xj) 

Xi.Xj^D 


wq^’IP‘ 


E exp (-^|LP - (jviTo 


X 6 ^E exp i-pUcore), (B22) 


where the <5 function enforces charge neutrality in the 
Green’s function representation. 


Appendix C: Polarization 


We consider the sum of each component of the electric 
field over the entire lattice in order to analyse the har¬ 
monic mode. The sum of the x /j/-component is split into 
separate sums over all cc/y-components that enter a par¬ 
ticular strip of plaquettes of width a that wrap around 
the torus in the y/x direction. Each component of the 
harmonic mode E^/y is then expressed in terms of the 
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charge enclosed along each of the strips of plaquettes: 


energy of the system given by 


xGD 

L—2a L—a 

=a {x + a) Y 

a:—0 y—0 

+ LaY {l- {^^y) 

y^O 
L — a 

+ LaYE.{l,y) (C2) 

y=0 

= -LYxYpi^) + EaY^^{^^y)^ (C3) 

x—a y—a y—CL 


(Cl) 


Exix+^,y^ -,y 


which follows from applying Gauss’ law to each strip of 
plaquettes that wrap around the torus in the y direction. 
The same argument holds for the y component, hence, 
the harmonic mode is given by 

E = --P + -^Wo, (C4) 

Co Leo 

where P := X)x 6D ^/^(^)/-^ i® origin-dependent 
polarization vector of the system and wq^x '■= 
£oaJ2y=a Exio,/‘2',y)/q is the x component of the origin- 
dependent winding field, with the y component defined 
analogously. Here, P and Wg are measured from a spe¬ 
cific origin. Note that the above applies to systems com¬ 
posed of either single- or multi-valued charges. 

We have thus shown that E, which is origin- 
independent, is given by the sum of two origin-dependent 
terms. One of these is attributed to the polarization 
of the system, while the other describes the winding of 
charges around the torus given that the polarization is 
measured with respect to the chosen origin. 

Restricting our attention to the gas of elementary 
charges, we now devise an origin-independent measure of 
the topological sector of the system. First, we note that 
adding lo windings to either component of the harmonic 
mode E corresponds to 

Ex/y ^ Ex/y + (C5) 

Leo 


and that this results in a change in the grand-canonical 


A[/ = ^oj ( -f 2E. 
Leo 


^xjy 


(C6) 


Hence, given an arbitrary charge distribution, the lowest- 
energy harmonic mode that describes the charge distri¬ 
bution is an element of the set in Eq. ( [To| . We therefore 
define a convention in which the harmonic mode is given 
by Eq. Q , where the polarization component of the har¬ 
monic mode is an element of the set in Eq. (101 and the 


winding component of the harmonic mode is given by Eq. 

Appendix D: Simulation details 


The system was simulated using the MR algorithm on 
an Lx L lattice of lattice spacing a = 1. One charge-hop 
sweep corresponded to picking a charge site at random, 
picking the x oi y direction at random, then proposing a 
charge hop in the positive or negative direction (at ran¬ 
dom) , repeating this 2N times (replacing each site / field 
bond after each proposal). One auxiliary-field sweep cor¬ 
responded to picking a charge site at random and propos¬ 
ing a field rotation around the site, repeating this N 
times. One global sweep corresponded to proposing a 
winding update in the positive or negative (at random) 
X and y directions. For all simulations, we performed 
five auxiliary-field sweeps per charge-hop sweep, and, for 
those simulations that also employed the global update, 
we performed one global update per charge-hop sweep. 
One charge-hop sweep corresponds to one Monte Carlo 
time step. 

The data sets in Sect ions [HT] and |TV] were averaged over 
multiple runs of 10® charge-hop sweeps per lattice site. 
The data set in Fig. was averaged over 608 and 446 
runs between T = 1.15 and 1.45 with the global update 
off and on, respectively, over 384 runs between T = 1.5 
and 1.6, and over 256 runs between T = 1.65 and 1.75. 

The L = 8 data set in Fig. |^was averaged over 128 
(T = 0.1-1.1), 256 (T = 1.15-1.39;T = 1.41-1.44;T = 
1.46 - 1.49), 768 (T = 1.4; T = 1.45; T = 1.5 - 1.75), 
and 256 (T = 1.8 — 2.5) runs; the L = 16 data set was 
averaged over 128 (T = 0.1 —1.1) and 256 (T = 1.15—2.5) 

runs; the L = 32 data set was averaged over 128 (T = 

0.1 - 1.1), 256 (T = 1.15 - 2.0), and 128 (T = 2.0 - 2.5) 

runs; the L = 64 data set was averaged over 128 (T = 

0.1 - 1.1), 448 (T = 1.15 - 1.45), 384 (T = 1.5 - 1.6), 
256 (T = 1.65 - 2.0), and 128 (T = 2.05 - 2.5) runs. 

We also simulated the L = 10, L = 20, and L = 40 
systems over small temperature ranges to calculate ad¬ 
ditional crossover points for Fig. all data sets were 
averaged over 512 runs. 
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